#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************        趋势分析讲解      *************
#*************                          *************

# 最初起源于剂量反应（Dose-Response）关联的分析
# 主要解决的问题是标准的分类分析没有有效利用类别内信息
# 回归方法通常用于测试趋势
# 通常列出每个类别对结果的影响强度(即比值比)与参考水平的 p 值以及“趋势”。


#### 0.准备好环境 ####
library(gtsummary)
library(survey)
library(haven)
library(tableone)
library(rms) # RCS-Regression Modeling Strategies
library(MatchIt) # PSM
library(plyr)
library(tidyverse)

setwd('G:/BaiduNetdiskDownload/NHANES')

# 运行一下读取 mort 数据的函数
NHSMortRead <- function(file.path){
  srvyin <- file.path   # full .DAT name here
  srvyout <- "data" # shorthand dataset name here
  
  # Example syntax:
  #srvyin <- paste("NHANES_1999_2000_MORT_2019_PUBLIC.dat")   
  #srvyout <- "NHANES_1999_2000"      
  
  
  # read in the fixed-width format ASCII file
  dsn <- read_fwf(file=srvyin,
                  col_types = "iiiiiiii",
                  fwf_cols(seqn = c(1,6),
                           eligstat = c(15,15),
                           mortstat = c(16,16),
                           ucod_leading = c(17,19),
                           diabetes = c(20,20),
                           hyperten = c(21,21),
                           permth_int = c(43,45),
                           permth_exm = c(46,48)
                  ),
                  na = c("", ".")
  )
  
  # NOTE:   SEQN is the unique ID for NHANES.
  
  # Structure and contents of data
  str(dsn)
  
  
  # Variable frequencies
  
  #ELIGSTAT: Eligibility Status for Mortality Follow-up
  table(dsn$eligstat)
  #1 = "Eligible"
  #2 = "Under age 18, not available for public release"
  #3 = "Ineligible"
  
  #MORTSTAT: Final Mortality Status
  table(dsn$mortstat, useNA="ifany")
  # 0 = Assumed alive
  # 1 = Assumed deceased
  # <NA> = Ineligible or under age 18
  
  #UCOD_LEADING: Underlying Cause of Death: Recode
  table(dsn$ucod_leading, useNA="ifany")
  # 1 = Diseases of heart (I00-I09, I11, I13, I20-I51)
  # 2 = Malignant neoplasms (C00-C97)
  # 3 = Chronic lower respiratory diseases (J40-J47)
  # 4 = Accidents (unintentional injuries) (V01-X59, Y85-Y86)
  # 5 = Cerebrovascular diseases (I60-I69)
  # 6 = Alzheimer's disease (G30)
  # 7 = Diabetes mellitus (E10-E14)
  # 8 = Influenza and pneumonia (J09-J18)
  # 9 = Nephritis, nephrotic syndrome and nephrosis (N00-N07, N17-N19, N25-N27)
  # 10 = All other causes (residual)
  # <NA> = Ineligible, under age 18, assumed alive, or no cause of death data available
  
  #DIABETES: Diabetes Flag from Multiple Cause of Death (MCOD)
  table(dsn$diabetes, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  #HYPERTEN: Hypertension Flag from Multiple Cause of Death (MCOD)
  table(dsn$hyperten, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  # Re-name the dataset, DSN, to the short survey name then remove other R objects
  assign(paste0(srvyout), dsn)
  rm(dsn, srvyin, srvyout)
  # View(data)
  return(data)
}

##### 1.1 提取数据模块 ##### 
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)
dim(demo.i) # 9971

# 提取生存数据 ：https://www.cdc.gov/nchs/data-linkage/mortality-public.htm

# 连接 National Death Index (NDI) ，补充 NHANES 参与人员的死亡数据
# 得到公开的 linked mortality files (LMF) 提供从调查参与日期到2019年12月31日的死亡率随访数据。
# 公开的数据中仅有随访时间或潜在的死亡原因
# 生存地址原始下载链接：https://ftp.cdc.gov/pub/Health_Statistics/NCHS/datalinkage/linked_mortality/
# 可以下载不同周期的生存数据，下载后进行读取

# setwd("C:/PUBLIC USE DATA")
mort.data_2015_2016 <- NHSMortRead('NHANES_2015_2016_MORT_2019_PUBLIC.dat')
# eligstat
# A value of 1 for ELIGSTAT indicates that the survey participant was eligible for the mortality linkage, 
# a value of 2 indicates the survey participant was under age 18 and not eligible for public release, 
# a value of 3 indicates the survey participant was not linkage eligible due to having insufficient identifying data to conduct data linkage.

# ucod_leading: 10种 潜在死亡原因，001-010，分别对应的原因看课程附件材料
# mortstat: 1生存状态
# permth_int: 从interview到随访时间点，间隔的月份数
# permth_exm: 从MEC-移动检查车到随访时间点，间隔的月份数

dim(mort.data_2015_2016) # 9971
# 转变为大写的列名，和 demo 进行拼接
colnames(mort.data_2015_2016) <- toupper(colnames(mort.data_2015_2016))
# View(mort.data_2015_2016)


##### 1.2 合并数据提取变量  & 去掉 NA 的行 & 衍生变量 ##### 
analyze.sample.data <- demo.i[,c('SEQN',
                                 "RIDAGEYR", "RIAGENDR", "INDFMPIR", "DMDEDUC2",
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]

analyze.mort.data <- mort.data_2015_2016[,c('SEQN', 'ELIGSTAT', 
                                            "MORTSTAT", "PERMTH_INT")]
# 拼接 demo & mort 数据
analyze.sample.data.add.mort <- join_all(list(analyze.sample.data, analyze.mort.data))

dim(analyze.sample.data.add.mort) # 9971    
# View(analyze.sample.data.add.mort) 

# 一键去掉 NA 的行
analyze.sample.data.add.mort.drop.na <- drop_na(analyze.sample.data.add.mort) 
dim(analyze.sample.data.add.mort.drop.na) # 5071

##### 1.3 衍生需要进行 trend 检验的变量 #####

# 衍生 Age 的三类变量
Age <- analyze.sample.data.add.mort.drop.na$RIDAGEYR
analyze.sample.data.add.mort.drop.na$Age.factor <- factor(ifelse(Age < 40, '1:Young', 
                                                                 ifelse(Age < 60 & Age > 40, '2:Middle', '3:Aging')))
# 不加序号，原始的排序：Aging Middle Young，按照首字母

##### 1.4 生成复杂抽样 NHANES_design ##### 
NHANES_design <- svydesign(
  data = analyze.sample.data.add.mort.drop.na, 
  ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTINT2YR, 
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

##### 1.5 趋势分析 #####
fit <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ Age.factor + RIAGENDR, 
                design = NHANES_design)
fit
tbl_regression(fit)

fit.category.trend <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ ordered(Age.factor) + RIAGENDR, 
                               design = NHANES_design)
fit.category.trend
a<-tbl_regression(fit.category.trend)%>%add_global_p()

show_header_names(a)
# .L 代表线性趋势
# .Q 代表二次趋势
